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Abstract 

In this paper, we extend the unsplit staggered mesh scheme (USM) for 2D magnetohydrodynamics (MHD) [D. 
Lee, A. Deane, An Unsplit Staggered Mesh Scheme for Multidimensional Magnetohydrodynamics, J. Comput. Phys. 
228 (2009) 952-975] to a full 3D MHD scheme. The scheme is a finite-volume Godunov method consisting of a 
constrained transport (CT) method and an efficient and accurate single-step, directionally unsplit multidimensional 
data reconstruction-evolution algorithm, which extends Colella's original 2D comer transport upwind (CTU) method 
[P. Colella, Multidimensional Upwind Methods for Hyperbolic Conservation Laws, J. Comput. Phys. 87 (1990) 446- 
466]. We present two types of data reconstruction-evolution algorithms for 3D: (1) a reduced CTU scheme and (2) a 
full CTU scheme. The reduced 3D CTU scheme is a variant of a simple 3D extension of Collela's 2D CTU method 
and is considered as a direct extension from the 2D USM scheme. The full 3D CTU scheme is our primary 3D solver 
which includes all multidimensional cross-derivative terms for stability. The latter method is logically analogous to 
the 3D unsplit CTU method by Saltzman [J. Saltzman, An unsplit 3D upwind method for hyperbolic conservation 
laws, J. Comput. Phys. 1 15 (1994) 153-168]. The major novelties in our algorithms are twofold. First, we extend the 
reduced CTU scheme to the full CTU scheme which is able to run with CFL numbers close to unity. Both methods 
utilize the transverse update technique developed in the 2D USM algorithm to account for transverse fluxes without 
solving intermediate Riemann problems, which in turn gives cost-effective 3D methods by reducing the total number 
of Riemann solves. The proposed algorithms are simple and efficient especially when including multidimensional 
MHD terms that maintain in-plane magnetic field dynamics. Second, we introduce a new CT scheme that makes use 
of proper upwind information in taking averages of electric fields. Our 3D USM schemes can be easily combined with 
various reconstruction methods (e.g., first-order Godunov, second-order MUSCL-Hancock, third-order PPM and fifth- 
order WENO), and a wide choice of ID based Riemann solvers (e.g., local Lax-Friedrichs, HLLE, HLLC, HLLD, 
and Roe). The 3D USM-MHD solver is available in the University of Chicago Flash Center's official FLASH release. 

Key words: MHD; Magnetohydrodynamics; Constrained Transport; Corner Transport Upwind; Unsplit Scheme; Staggered Mesh; High-Order 
Godunov Method; Large CFL Number 



1. Introduction 



Many astrophysical applications involve the study of magnetized flows generating shock waves. Model- 
ing such flows requires numerical solution of the equations of magnetohydrodynamics (MHD) that couple 
the magnetic field to the gas hydrodynamics using Maxwell's equations. A valid computer model needs to 
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capture accurately the nonlinear shock propagation in the magnetized flows without sacrificing computa- 
tional efficiency and stability. 

Obviously, with suitable assumptions about flow symmetries, a simple approach to obtain a computation- 
ally efficient model is to solve a reduced system in ID or 2D instead of 3D. However, a limitation of such 
reduced systems is that they cannot be used to understand complicated nonlinear physics occurring only in 
the full 3D situation. Although solving the reduced system can illustrate interesting characteristic features 
(e.g., the inverse energy cascade in 2D turbulence [39]), it is essential to use 3D simulations term in order 
to understand the full nonlinear nature of MHD phenomena (e.g., the energy cascade from large scales to 
small scales in 3D turbulence) . 

There are two approaches in modeling multidimensional (i.e., 2D and 3D) algorithms for gas hydrody- 
namics and MHD in terms of spatial integration methods: split and unsplit. The directionally split method 
has the advantage of extending a ID algorithm to higher dimensions, simply by conducting directional 
sweeps along additional dimensions, in which each sweep solves ID sub-system. Thus, the Courant-Friedrichs- 
Lewy (CFL) numerical stability constraint of the split schemes in multi-dimensions is the same as the ID 
constraint, which is to say CFL< 1.0. Despite their simplicity and robustness, however, a number of re- 
cent studies have revealed numerical problems in the split formulations of multidimensional MHD and gas 
hydrodynamics (e.g., loss of expected flow symmetries [2,41], failure to preserve in-plane magnetic field 
evolution [30,40], numerical artifacts due to a failure to compute proper strain rates on a grid scale [3]). 

For MHD the use of an unsplit formulation is more critical than for hydrodynamics. This is because 
the split formulations fail to evolve the normal (in the sweep direction) magnetic field [18, 30, 31, 59]. For 
2D MHD, Gardiner and Stone [30] identified the importance of such multidimensional consideration in 
their unspUt MHD scheme based on the comer transport upwind (CTU) [16] and the constrained transport 
(CT) [26] methods. Later, the authors proposed a 3D unsplit version of an unsplit MHD scheme in [31], in 
which the extension of the multidimensional MHD terms from their 2D algorithm to 3D is accomplished at 
the cost of considerable algorithmic complexity and a reduced stability limit (CFL < 0.5) in their 6-solve 
CTU-i-CT algorithm. It is known in a CTU-type 3D unsplit formulation that the full CFL stability limit (i.e., 
CFL number < 1 .0) can be recovered by accounting for intermediate Riemann problems fully, requiring 12 
Riemann solves per zone per time step [55]. In general, the calculations associated with the Riemann solves 
are computationally expensive. Gardiner and Stone [31] considered two alternative options, an expensive 
12-Riemann solve yielding the full CFL limit and a reduced 6-Riemann solve with a more constraining CFL 
condition (CFL number < 0.5). They found that the two approaches are similar in terms of computational 
cost and there is no significant difference in performance between them. The 6-solve scheme is chosen to be 
their primary 3D integrator because of its relatively low complexity in incorporating the multidimensional 
MHD terms. 

The CTU formulation has an advantage in its compact design of one-step temporal update which is 
well-suited for multidimensional problems. However, it is limited to second-order. There has been much 
progress in other types of temporal update strategies that are higher than second-order accurate, taking a 
different path from CTU. Early attempts have utilized a Runge-Kutta (RK) based temporal update formu- 
lation coupled with spatially high-order reconstruction schemes in the finite-difference framework [6, 13, 
35,38,44,57,58,60]. Such RK-based high-order schemes have been also developed in the finite-volume 
framework [8, 22, 23, 37, 69] which has superior properties to that of finite-difference for resolving com- 
pressible flows on both uniform and AMR grids. The high-order RK temporal update strategies rely on 
multi-stage updates which add to the computational cost. Therefore it is desirable to retain a CTU-like 
one-step formulation, while retaining higher than second order accuracy. Recent work has been found 
to provide such efficiency using a new formulation so-called the Arbitrary Derivative Riemann Problem 
(ADER), see [9, 12, 24, 25, 61, 62, 64]. For solving multidimensional conservation laws, there has been an- 
other line of progress that tries to build genuinely multidimensional Riemann solvers for hydrodynamics 
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[1, 15,27,28,32,67]. Recently, a family of two-dimensional HLL-type Riemann solvers, HLLE [10] and 
HLLC [1 1], have been introduced and generalized by Balsara for both hydrodynamics and MHD. As shown 
in his work the multidimensional Riemann solvers are genuinely derived for 2D. A major improvement in 
MHD flows is that they inherently provide proper amount of numerical dissipation that is necessary to prop- 
agate magnetic fields in a stable manner. Alternatively, ID Riemann solver formulations such as [30,45] 
need to add extra dissipation for a stable upwinding. The use of multidimensional Riemann solvers is also 
shown to capture isotropic wave propagations better than the usual ID approach. Furthermore, both types 
of solvers have been extended to 3D using a one-step predictor-corrector formulation. 

The above mentioned strategies using high-order schemes and genuinely multidimensional Riemann 
solvers, provide improved solution accuracy and stability over CTU-CT formulations. In this paper, how- 
ever, we are primarily interested in constructing a scheme that can be built on the ID Riemann solver frame- 
work in line with a CTU-type method. The latter is (arguably) most widely used in many Godunov-type 
modern codes [16,30,31,40,43,45-47,55]. This design also benefits us in extending our 2D USM-MHD 
algorithm [40] to 3D without any modifications of the Riemann solvers. This paper describes an approach 
that provides (i) an algorithmic extension from 2D to 3D of the USM scheme of Lee and Deane [40], and 

(ii) the full CFL stability bound in 3D without the expense of 12 Riemann solves per cell per time step, and 

(iii) a new upwind biased electric fields construction scheme for CT. We show that the present USM scheme 
achieves a numerically efficient and consistent MHD algorithm in 3D without introducing a greater amount 
of additional complexity, while maintaining the full CFL stability range. 

The paper is organized as follows: Section 2 describes our new 3D unsplit, single-step data reconstruction- 
evolution USM algorithm which consists of two stages, i.e., normal predictor and transverse corrector. Sec- 
tion 2 is subdivided into several subsections. We begin in Section 2. 1 our discussion of the 3D USM scheme 
by considering the governing equations of MHD and their linearized form. The second-order MUSCL- 
Hancock approach for calculating the normal predictor is described in Section 2.2. We introduce in Section 
2.3 our two 3D CTU schemes to compute the transverse correctors, which are efficient and essential for 
obtaining the full CFL stability range. In the subsections therein, we construct Riemann states at cell inter- 
faces, focusing on our new transverse correctors that do not require the solution of any Riemann problem. 
The Riemann state calculations are completed by evolving the normal magnetic fields by a half time step, 
about which we describe in Section 2.4. The final update of the cell-centered conservative variables is 
shown in Section 2.5, followed by a new 3D upwind-biased CT update of magnetic fields in Section 2.6. 
We summarize our step-by-step, point-to-point 3D CTU schemes in Section 3. In Section 4 we present 
numerical results of various test problems that demonstrate the qualitative and quantitative performance of 
our schemes. We conclude the paper with a discussion in Section 5. 

2. The three-dimensional USM scheme for MHD 

2.1. MHD Equations 
We consider solving the equations of MHD in conservation form 




(1) 



(2) 



(3) 



(4) 
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The conservative variables include the plasma mass density p, momenta pu, magnetic fields B, and total 



energy density E. The rest are the thermal pressure p = [y— \){E 



ip[/2 



Bp), the magnetic pressure 



Bp = {B^ + B^.+B,) /I, and the sum of the two is the total pressure p,o, = p+Bp. The ratio of specific heats 



is denoted with y as usual. The solenoidal constraint V • B 
We write the above equations in a matrix form in 3D 



is implied in the induction equation. 
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(5) 



where U contains the eight MHD conservative variables, and F, G, and H represent the corresponding 
conservative fluxes in x,y and z directions. It is often convenient to cast the conservative form of Equation 
(5) into a quasi-linearized representation in terms of primitive variables, V =(p , M , V, w, Bj , By , Bj- , p ) , in 
order to discretize the coupled system of MHD equations (l)-(4), 

av 



av av av 
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The coefficient matrices A^, A,,, and A- are given by 
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(6) 



(7) 



(8) 



(9) 



For exposition purposes in this paper, we illustrate our calculations using a spatially second-order MUSCL- 
Hancock (MH) piecewise-linear method (PLM) for the normal predictor. Other normal predictor algorithms 
(e.g., piecewise parabolic method (PPM [17]), essentially non-oscillatory (ENO [34]), weighted essentially 
non-oscillatory (WENO [38]), etc.) can be adopted as well to give different degrees of solution accuracy 
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in our algorithm. In fact, we have implemented various reconstruction schemes of MH, PPM and 5th order 
WENO in FLASH, and they are available in the official FLASH distribution [21,29,40,41]. 
This brings us to write the system (6) to obtain second-order accurate discretizations at cell faces, 

V- tlS^w = n,t + ^ [±I - K,,, - ^^y^Vj, - ^A,Af V;:,„ (10) 

V- = V;:,, - ^A.A-V;: + i [±I - ^A„] Af V;: - ^A,A-V«,„ (11) 

v" ti:?. = y'lj, - ^A.A^^'^M " ^y^y^rnj,^ +^ [±i - ^aja-^v;;,,, (12) 

where the plus and minus signs correspond to directions of N,E,S,W,T and B respectively in a natural way, 
see Figure 1 . Each matrix represents the coefficient matrix in the <f -direction evaluated at V" y ^. The 
undivided difference operators in each J-direction are denoted as A'J'^ and A^'^ and they are suitably chosen 
slope vectors of V"^^ in each cell {i,j,k) using TVD and upwind slope limiters, respectively. 





V"+l/2 








V"+l/2 








V"+l/2 


V"+'/2 
^ i+l.j.W 




V«+l/2 





Fig. 1. The boundary extrapolated values on a 2D cell geometry. Our subscriptions N,S,E,W,T,B represent respectively north, south, east, west, 
top and bottom that are based on a reference point at the local cell center node {i,j,k). 

2.2. Normal Predictor 

The first stage is to calculate the normal predictor states, including all the required multidimensional 
MHD terms (the MHD terms hearafter) [40] satisfying the solenoidal constraint V • B = 0. We begin our 
discussion with the evolution of the normal field, B^, which is treated separately from the other primitive 
variables. For instance, when N = x,'we can define 





V, 




A.v As, 
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and A V = 












(13) 



Here Vjc is a 7 x 1 vector excluding Bj^, A;r is a 7 x 7 matrix omitting both the fifth row and column in the 
original matrix A^ in Equation (7), and A^^ is a 7 x 1 vector. 
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U, , , , —V, — w, — /cu • B 

P P P 



(14) 



Note that the hat {") notation denotes the reduced system (i.e., the one corresponding to the usual ID MHD 
equations) and the bar (-) notation indicates the re-assembled full system. Similarly for the other directions, 
we have 
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The term Ag^ for each will be our representation of the corresponding MHD term in this paper. 

The first step of MH extrapolates V" ^ ^. to construct the six multidimensional Riemann states V" j^^^ sewtb 
at cell interfaces to achieve second-order accuracy by using a total variation diminishing (TVD) slope lim- 
iter * . Although the slope limiter can be applied to either primitive or characteristic variables, we prefer 
the latter since it is less prone to generating spurious oscillations as noted in the literature [63, 66]. We do 
not apply any limiting to Bj^, allowing the continuity of the normal field at cell faces (e.g., see discussion 
in [40]). To simplify our discussion, we focus on the A:-direction in Equation (10). The others in Equations 
(1 1)-(12) can be computed in the similar way. We consider the first two terms in (10) that are related to the 
normal predictor 
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(17) 



where A^^Vj; • 



and AB'l - 



... The notations B,j and denote 



x,/+l/2,,/,yt " ""x.i-l/lj.k- "-tcttiv^iio dnu Dd 
cell-centered and cell face-centered magnetic fields respectively, with d = x,y,z- In CT, AB" , is constructed 
such that the numerical divergence is zero using the cell face-centered magnetic fields. In other words, AB" ■, 



AB" I and AB"^, are chosen such that 

AB". AB". AB" 
Ak Ay Az. 



(18) 



where we analogously define AB"^ and AB"^. Solving a system in relation (17) is equivalent to considering 
two sub-systems 



The second relation in (19) is nothing but 

/n \n+l/2,|| _ nil I ^ AD" _ ^« 

\^^)iJ,k,E,W — ^x,ij,k^^^x,i — '^x,i±l/2,j,k^ 



(19) 



(20) 



* For instance, limiters such as minmod, van Leer's, monotonized central (MC), or a hybrid coinbination of them on different wave structures [7] 
can be used. 
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because we use a simple arithmetic averaging to obtain the cell-centered magnetic field using the divergence- 
free fields at cell interface centers, 



^",i,j,k — 2 {^'xJ+l/2,j,k~^K,i-l/2,i,k J ' 

Applying the characteristic tracing method in ;c-normal direction in (19) yields 
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A suitable TVD slope Umiter along the ;<;-normal direction is used in the undivided slope operator on each 
characteristic variable a 



A';''aljj^ = TVD_Limiter 



i+l,j,k' 
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(24) 



Here X'"- ^ j.,r'", j k^K\ j k represent respectively the eigenvalue, right and left eigenvectors of Av, calculated 
at the corresponding cell center {ij,k) in the .r-direction at time step n. This completes the first part of our 
description on a single-step, data reconstruction-evolution algorithm in the ;c-normal direction. 



2.3. Transverse Corrector in USM 



2.3.1. Review of Computing Transverse Flux Gradients using Characteristic Tracing 

The transverse corrector adds the gradients of transverse fluxes to the normal predictors. This transverse 
corrector step plays a crucial role for stability in CTU. Generally speaking, the degree of accuracy is affected 
by the normal predictor, whereas numerical stability is strongly determined by the transverse corrector [14]. 

In [40], Lee and Deane noted that the transverse flux gradients, responsible for the cross-derivative terms 
in CTU, which assure stability for flows advecting along diagonal comer directions, can be replaced by a 
simpler approach that is based on characteristic tracing alone. This removes the need to solve the intermedi- 
ate Riemann problems. As a result, this approach requires only two Riemann solutions in 2D (not counting 
the extra two Riemann solves to update the divergence-free magnetic fields by CT), while preserving the 
full stability of the CTU scheme. We review a pointwise description of the transverse corrector in USM for 
a moment. Consider the j-transverse flux gradient (i.e., the third term in (10)) which supplies the corrector 
term to the ;c-normal predictor states. For any left (V/ y-,/.) and right (V/,y+ ^.) states at cell {i,j,k) along 
3'-direction, the jump conditions across the individual m-th wave gives 

"lo-l ,^ ^ ^ 

^y,i,j,kyy,i,j-,k+ ^ ^y]ij,k^'^,ij,k^y''^lj,k = ^y,i,j,k^yJ,j+,k ^ ^ ^'y,ij,k^y,ij,k^y''^'ij,k- (^5) 
m= 1 m=mo 

Now recall that the property of conservation [42, 63] across discontinuities of the Roe matrix A. It states 
that the Roe matrix ensures conservation across a discontinuity between the left (V/) and right (V^) states, 
given by Flux(V^) — Flux(V/) = A(Vr — V/). Applying this relation to Ayjj^k, ^yjj-.k and V,,,,-.y+i:, and 
from (25), we obtain 
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Gij+i/2 - G/,/-i/2 = Ayjj^k{\yjj+^k - Vv./j-i) = Xyijf,r^iji^A"Pa"jj^. (26) 
The upwind slope limiter A'C is applied to each characteristic variable a" j as 
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In Equation (26) we see that the sum over all wave contributions gives an effective upwinding of trans- 
verse flux gradients in j-direction. The advantage in this approach is that there is no need to solve the 
intermediate Riemann problems to add the transverse flux gradient correction terms to the spatially re- 
constructed, temporally evolved, normal predictor states in order to gain the upwind stability in the CTU 
formulation. Because we rely on using the eigensystem in Equation (26), one might suspect that this char- 
acteristic tracing approach could be as expensive as directly solving the associated Riemann problems at 
each interface, followed by taking the gradient of the computed transverse fluxes. However, this is not the 
case because we reuse the y- (or x-) directional eigensystems that were already calculated in the normal 
predictor step in the y- (or x-) direction. Thus there is no need to compute any additional eigenstructure for 
each transverse direction, which makes our scheme much more computationally efficient than the standard 
CTU method. A Fortran-like pseudo code illustrating the algorithm is as follows: 

do j=jmin, jmax 

do i=imin, imax 

! Compute normal predictor in x-direction, and 

! store x-directional normal predictor states & eigensystems in arrays 

call dataReconstructNormalDirection (x_dir, x_normalPredictStates, sigmaSum_x) 

! Compute normal predictor in y-direction, and 

! store y-directional normal predictor states & eigensystems in arrays 
call dataReconstructNormalDirection (y_dir, y_normalPredictStates, sigmaSum_y) 

! Transverse Correction to the x-normal predictor 

x_normalPredictStates = x_normalPredictStates - . 5*dt/dy*sigmaSum_y 
! Transverse Correction to the y-normal predictor 

y_normalPredictStates = y_normalPredictStates - . 5*dt/dx*sigmaSum_x 



In the above, the terms sigmaSum_x and sigmaSum_y represent the summation of all wave contributions in 
the X- and j-directions, respectively, given in Equation (26). The rest of the terms are self-explanatory. 

Our approach to approximate the transverse flux gradients, solely using the characteristic tracing, greatly 
simplifies the overall unsplit CTU algorithm by reducing the number of required Riemann solves. In gas 
hydrodynamics, the proposed algorithm requires a total of three Riemann solves to update the solution from 
n to « + 1 without compromising solution stability and accuracy. It will be shown for MHD in Section 2.4 
that three additional Riemann problems (yielding a total of six) are required to update the divergence-free, 
cell face-centered magnetic fields using the CT method. Another advantage in our approach, especially for 
MHD, is the relatively simple handling of multidimensional MHD terms. This is because our method of 
adding transverse flux gradients provides a single-step, directionally unsplit data reconstruction-evolution 
algorithm to calculate Riemann states at cell interfaces. It is therefore much simpler to enforce the balance 



end do 



end do 
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between flux gradients in all three directions associated with the MHD terms. As noted in [31], complica- 
tions arise in the standard full 12-solve CTU scheme, in which the MHD term balance seems to be hard to 
achieve in a series of partial transverse flux gradient updates based on dimensional splitting. 

2.3.2. Reduced 3D CTU Scheme in USM: Interface State Update from n to n + \ /2 Time Step 

Our first simple algorithm using the transverse corrector technique in the previous section is analogous 
to the 6-solve CTU in [31]. This approach can be viewed as a straightforward 3D extension of the 2D 
CTU scheme [16], omitting all the third-order cross-derivative terms such as /d^d^d,, while including 
the second-order cross-derivative terms that are provided in the 2D CTU method. The resulting Riemann 
state calculations account for flow information along the edge directions, but do not fully account for flow 
information along the diagonal comer directions, yielding a formal stability limit of CFL number less than 
0.5 ^ . This simple approach, referred to as the reduced 3D CTU scheme, can be directly extended from the 
2D CTU [40] by adding the third additional transverse flux correction in z- That is, the .x-normal predictors in 
Equations (22)-(23) are further corrected by including the transverse flux contributions from j,z-directions 
using the characteristic tracing approach described in the previous section, see also [40]. For instance, in 
Equation (10) the transverse corrector step can be updated, first by accounting for the y-trans verse flux 
correction, 

y«+l/2,y _y,+l/2,|| _ ^ A CV« U^'^V" QS^ 
^ i,j,k,E.W ~ ^ i,j,k.E.W 2Ay - 

followed by the z-transverse flux correction. 

In these transverse corrector steps, it is important to use the upwind biased slope limiter instead of any form 
of TVD limiters as reviewed in Section 2.3.1. Note that in the original 2D CTU scheme by Colella [16], 
using the upwind flux gradients in the transverse directions is the key mechanism that guarantees the full 
CFL stability bound. We establish the same upwind couplings by means of using the upwind slope limiter 
for our transverse corrector. Using a TVD slope limiter instead would result in a reduced stability limit for 
our algorithm (and we avoid using it). The two transverse correction terms in (28) and (29) are calculated 
as in Section 2.3.1, completing our description of the reduced 3D CTU scheme. 



2.3.3. Full 3D CTU Scheme in USM: Interface State Update from n to n + 1/2 Time Step 

To establish the full stability limit (CFL number less than 1 in 3D) as featured in the 12-solve CTU 
scheme of Saltzman [55], we need one more step to couple diagonally moving flow effects. This situation 
occurs when the conservative quantities are advected across the corners diagonally with components of 
the local velocity fields (m,v,w) being of comparable orders of magnitude. In USM, these couplings can 
be added to the interface states by performing intermediate state calculations at ?i + |. They involve extra 
evaluations of the coefficient matrices and the undivided upwind differences in (28) and (29) at 

y«+i/3,z_y„ ^(\\" ,A"''v" , no'i 

and 

v^tf = v^,,. - ^i^^hkK'nj,- (31) 



' One can easily prove this stability bound numerically for a 3D scalar advection equation using a standard von Neuman Fourier analysis, assuming 
a single Fourier mode solution g'l j k = e'(y^+1-'^>'+?*'^) where i = \/— 1; I,J,K as the grid indices; and ^,T|, t, the wave numbers in j:,y,z-directions 
respectively. 
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More specifically, the transverse correctors in Equations (28) and (29) are replaced by 

^ i,j,k,E,W ~ ^ i,j,k,E.y/ 2Aj ''■''''^ ^ ''^'''^ ' 



and 



(32) 



(33) 



Here we make one important observation. Note that the additional re-evaluations of the matrices Ay and 
A-, at the « + ^ states V"!^^'' and V"|],^^'^ simply mean that the corresponding eigensystems for the charac- 
teristic tracing in the transverse directions need to be re- calculated, incurring the corresponding additional 
cost. Considering the full 3D interface state calculations in Equations (10)-(12), there are a total of six 
additional eigensystem evaluations required for the transverse correctors, which becomes as expensive as 
directly solving the corresponding Riemann problems, making our scheme expensive. Therefore an efficient 
alternative approach is required. Noticing 



— (A)" 



(34) 



V".. 



and using a Taylor expansion at V" ^ ^, we consider 



3G 




dG 


M an 


a^G 




3V 




av 




av2 


V" 



(35) 



Ignoring the error term in the matrix evaluations in Equation (35), we can replace respectively Ay(V - ^ 



«+1/3,<-;n 



and Aj,(V"|y^''*) with Ay(V"y^) and Ayiy^^.j^ in Equations (32)-(33). However, it is essential to retain 



d\ dy 



, AvandAi'/'V 



■«+l/3,V 

i-j-.k 



^ J 1 Az to couple the diagonal upwind 



corner transport. We proceed this as follows. Ignoring the 0{At) term and keeping the first-order approxi- 
mation in Equation (35) for the matrix evaluation, the transverse corrector in Equation (32) becomes 



Using our transverse corrector strategy, we get 



^n+l/2,y _^«+l/2,|| 
^ y,i,j,k,E,W ~ ^ 



where the upwinding slope applied to each characteristic variable a is given by 



fT)i+l/3,z ,7«+l/3,z 



.j+\,k 



Notice that the MHD term at « + 4 can be written as 



f^«+l/3,z fwj+l/3,c 



) if K,J.<^ 
i,-i.k) if K.,.k>^ 



ij,k 



(36) 



(37) 



(38) 



AB 



yj 



Ay ( b;,, - ^ KA^lj^tK'%^k + i^EXj^t^lk] ■ e«,. ) = Ay ( + 0{At] 



(39) 



where e^, is a unit vector in By direction for contraction and the hat notation implies the omission of the B^ 
components. However, in order to choose AB"'^j^^^'^ to enforce the numerical divergence to be zero always 
(see Equation (18)), we further drop the At error term and only take 
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■ nn+1/3,; _ . n„ _ I « _in 



(40) 



'y,,/ -v"y.7 "y^z+V^ "yj-i/S' 

where Z?" ^^j^j ^i"^ the cell face-centered, divergence-free magnetic fields in y-direction. 

The upwind differences in relation (38) are given by (assuming uniform spacing in each direction every 
where), 



^«+l/3,z_^n+l/3,z 
^ i,j,k 

- V" _ V" ^ 

-^i,j+X,k ^i,j,k 3^ 

-V" -V" - — 

— ^ i.j+\.k ^ ij.k 3^ 



12 ^z.i.i+l.k^'z,i,i+\-,k^'z'^'i,i+l.k ~ ^ ^[.L i.k^'z,i,j,k^T^l i,k 

1=1 



h=l 



+{AlABl),j+i^k-{AlAB",)ij, 



(41) 



and 



frn+l/3.z _^n+l/3,z 
^ij,k ^ij-l,k 

- V" - V" - — 

— V" — V" — 



7 7 
XI ^Z:iJ:k^Z;iJ,k^z^^'i,j,k ~ ^ \,i,i-\,k^ Z;i;i-\,k^z' ^'i,i-\-,k 



ij-i,k 



h=l 



1=1 



In the case of X"!- , /. > for all m, using (42), Equation (37) becomes 

At ' 



(42) 



Af 



^«+l/2,y _-(V«+l/2,|| ^' \ T,m „m im /^)7+l/j,; fx \n A R" 



^«+l/2,|| V T,m „m im ("(Vn \ ^ tx \n \T>n 

^y,Lj,k,E,w ^ L\iJ,k^yJjry,i,j,k-\^ij,k ^ ij-i,k) 2Ay^ 



(44) 



6AvA7 i ^ ^y,lj,k^y.ij,l^^ij,k ' L Ki,j,kKi,j,k^T^i.,j.,k " L 
\m=l h=i 1=1 



l,k 



+{AlAB:)„^k-iAlABl)i, 



i-i.k 



(45) 



Note that the terms in relation (44) are what we already established in the reduced 3D CTU scheme. The 
terms in relation (45) are new correction terms for the full 3D CTU scheme that need to be added to the 
reduced 3D CTU interface states in relation (44). 

Similarly, the final form of .x-interface states w Equation (33) is established by adding another 

correction term appearing in A"p\"^1^^'^ , giving the result (assuming X"\ ^ ^ > for all m) 



^)7+l/2,y 



At J- 



' iJ,k,E,W ~ ' y,i,i,k.E,W 2Az 



\ij,k^z,i,j,kh,i,j,k' ij,k- ^ij. 



k-l, 



m=\ 



At 
2Ai 



{^B.)h,k^l, 



+ 



At'- 



^ 1 



(a«,,as;),-,,,-(a«^,ab;),-,,,_i 



52 ^'z,iJ,k^,iJ,kKi,j,k ' H ^y,i,j,k''y,i,j,k^y'^^'i,j,k ~ ^ '^'y,iJ,k-l''^y,iJ,k-\^T^'i,j,k-l 



■h=l 



1 

L 

1=1 



(46) 



(47) 
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Likewise, the terms in relation (47) are the extra correction terms required for the full 3D CTU scheme. 
They must be added to the reduced CTU terms in relation (46). 

It is worth pointing out at this stage that all of the eigensystems in the full 3D CTU correction terms 
are readily available as they have been calculated and stored in the normal predictor step in each x^y^z- 
direction described in Section 2.2. In the normal predictor step, one can store not only the eigensystems, but 
also the two summations inside the square brackets in relations (45) and (47) (see also the simple pseudo 
code in Section 2.3.1). The only extra calculations imposing additional computational costs are therefore 
the upwind differencings inside the square brackets and the dot products in relations (45) and (47), which 
are computationally much more efficient compared to the calculation requirements in the 12-solve CTU 
scheme. This completes our description of the single-step, data reconstruction-evolution algorithm for all 
variables except the divergence-free normal magnetic fields at each cell face. The reconstructed interface 
states are second-order accurate in space and evolved to « + | time step at each interface. The next step is 
to advance the remaining normal fields at the cell faces, finalizing the Riemann state calculations. 

2.4. Advancing the Normal Fields from n to n + I /2 Time Step using CT 

In updating the normal fields to ?i + ^, it is important to meet two conditions. The first is a continuity 
restriction of the normal magnetic field across cell interfaces [7, 18,30,40,51]. The second is the divergence- 
free constraint of the normal fields on a computational grid. As a last step of our Riemann state calculations, 
we must evolve the normal field components at each cell boundary by a half time step, while satisfying the 
two conditions. We therefore follow the CT approach using the high-order Godunov fluxes that are solutions 
to a Riemann problem using the Riemann states e w t b described in Sections 2.3.2 and 2.3.3. Our 

approach here is the 3D extension of the 2D method in using the same approach as in [40]. Here we briefly 
describe the procedure only in jc-direction, which can be similarly applied to the other directions. We first 
solve Riemann problem at x interfaces as 

* i-\l2.j.k - yi-l.j.k.E^ ^iJ,k.W ) ' * i+\/2.j.k - [^iJ,k.E ' i+lJ,k,W ) ' ^^S) 

With these high-order Godunov fluxes at the half time step we evolve the normal fields by a half time step 
using the CT update 

, «+l/2 _ , „ 
'^x,i+l/2,j,k — '^x,i+l/2J,k 

__^/f*'"+V2 _p*.n+l/2 \__^/ _ F*'"+V2 ,p*,n+l/2 \ ..qs 

2A3;l '•■'+'/2'.'+V2,<.- ^z.i+l/2J-l/2.kj 2Az\ v./+l/2,/,/:+l/2^-^y,/+l/2,y,«.-l/2/' ^^^^ 

where the duality relationship between the electric fields and the high-order Godunov fluxes [5] is assumed 
in the expression. The electric fields ' ^ in (49) can be constructed based on the MEC method [40] 
that takes an arithmetic average of four Taylor series expansions of the fluxes given by (48) to obtain them 
(e.g., see Equation (53) in Section 2.6.1). The normal fields in (49) satisfy the divergence-free constraint as 
well as the continuity restriction across cell interfaces as they are direct solutions to numerical induction 
equations via the CT approach. 

Given these updated cell face-centered divergence-free fields, the Riemann states at ;ic-interfaces are up- 
dated as 

^i.JXE ■ - '^xj+l/2,jr ^UXW ■ «Bv - '^xJ-l/2.J.V <-5") 

where e^^ are unit vectors for the magnetic field components in .r-direction. 



* Note here that we use a consistent superscript (e.g., F* and £*) between the Godunov fluxes and the electric fields that are in the duality 
relationship. The superscript is used for the intermediate Riemann solutions in Section 2.4, whereas the superscript * (e.g., F* and £") is used for 
the final Riemann solutions in Section 2.5. 
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Now that the second-order accurate Riemann states V" j"}^^ sewtb ^''^ available, the second-order Go- 
dunov fluxes can be evaluated by solving the last set of Riemann problems at x-interfaces, 

'^i-l/2J,k — \ ^ i-\.jAE^ ^ iJ.k,W J ' '^i+l/2J,k — i,j,k,E^ ^ i+lJ,kW ) ^ ^^^> 

Note that the superscript * is used to represent the second-order Godunov fluxes that are the solutions of 
the Riemann problems. 

2.5. Cell-centered Solution Update from n to n + l Time Step 

The USM algorithm updates the cell-centered conservative variables to the next time step n+l using an 
unspht integrator, 

¥T«+1_IT" ^ f|7.*,"+l/2 _|-,*,«+l/2 1 _ rp*,«+l/2 _p*.«+l/2 \ _ fTi*,n+l/2 _ti*.«+1/2 \ 
^iJ,*-^;,M^\*/+l/2,M *'/-1/2.,M/ ^y\^Lj+\l2,k ^/„/-l/2,i/ ^ \"/,M+l/2 ",„M-1/2/ -(^^^^ 

In general, after this update, non-zero divergence magnetic fields are still present at cell centers, and they 
need to be corrected. In the next section we update the divergence-free cell face-centered magnetic fields 
from nto n + \ time step using the modified electric field construction (MEG) scheme [40]. The cell face- 
centered fields are averaged to correct the cell-centered magnetic fields at the n + l state. The choice of a 
time step Af for the full 3D CTU scheme is limited by the full CFL bound, with which we set our CFL 
number to be 0.95 for all numerical results presented in this paper, unless otherwise stated. 

2.6. Face-centered, Divergence-Free Fields Update via CTfrom n to n+\ Time Step using MEC 

2.6. 1 . The Standard Arithmetic Averaging Approach in MEC: standard-MEC 

In [40], the 2D version of the modified electric field construction (MEC) scheme was introduced. The 
method provides electric fields at cell corners using high-order Taylor expansions. Displaying the electric 
field in z-direction only, this standard-MEC algorithm gives 



pn+1/2 _ p*,"+l/2 I Ay '^^.v+i/2,j,t , Ay- ''"'fe;.,+ i/2,j,t , ni' Ai,3^ 

■^Z,i+l/2,/+l/2,<: - ■^2,;+l/2,,/,*:^ 2 ay 8 ^ ^K^J ) 

;ic'+,n+l/2 -^2p*./I+l/2 

pn+1/2 _ p*,n+l/2 Ay '^^;,,+i/2j+i,t , Ay^ ''"-^;.i+i/2j+i.t , 

^z,i+l/2J+l/2.k~-^z,i+l/2J+l,k 2 dy "^8 ay^ -T uyisy 



;)r*-"+i/2 ^ ;i2f«,"+i/2 (53) 



F«+V2 _ p*,«+l/2 I Ax '^^.v.j+i/2,t I Ax' ^:.ij+i/"2,t I n( \^^\ 

^z,i+l/2,j+\/2,k- -^z.i.J+l/2,k^ 2 dx 8 dx^ -\- U{nx ) , 

pn+l/2 _ p*,«+l/2 Ax ''^:,/+l.,;+l/2,t I Ax' ^:..+ lJ+l/2.t I n( \y^\ 

.^z,i+l/2J+\/2.k~ -^zj+lj+l/2,k 2 dx "'"8 dx^- -\-U{lXX). 

The duality relationship [5] has been assumed for those electric fields at cell face centers about which the 
Taylor series are expanded. The standard-MEC algorithm proceeds to take a simple arithmetic average of 
these four Taylor expansions of each field in Equation (53) to get an averaged electric field E"'^^^^^ ;+i/2 k- 

2.6.2. Upwind Biased Averaging in MEC: upwind-MEC 

The standard CT approach of taking the arithmetic average of the four electric fields is simple enough 
to work well in local smooth regions. The simplest form of this averaging approach was first suggested 
by Balsara and Spicer [5] using ID based Riemann solvers. This idea seems a very natural choice if one 
considers the grid locations of the electric fields. However, the authors understood that in truly multidi- 
mensional flows where there is a directional bias, the simple arithmetic averaging scheme may need to 
be corrected and it would be better to incorporate upwind information. Recently, a general resolution on 
such issue with multidimensional upwinding has become available by the subsequent efforts to build the 
genuinely multidimensional HLL-type Riemann solvers by Balsara [10, 11]. 
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On the other hand, within the ID Riemann based CTU approach, Gardiner and Stone [30] identified the 
shortcomings of the simple arithmetic averaging method and developed a systematic construction of CT 
algorithms that are consistent for a plane-parallel, grid- aligned flow. They recovered the necessary amount 
of numerical viscosity that stabilizes their underlying CT integration algorithm. The CT methods proposed 
therein readily satisfy planar symmetry for d/dx = or d/dy = 0, showing that their algorithms recover 
the associated one-dimensional solution for the underlying integration algorithm. A similar approach of 
increasing dissipation is also found in [45]. 

Although the approach by Gardiner and Stone provides consistency for plane-parallel, grid-aligned flows, 
the method does not take into account multidimensional effects where the flow has one specific directional 
bias without assuming d/dx = or d/dy = 0. To illustrate this, we consider the weakly magnetized field 
loop advection problem [30, 40] where the loop is advected by a dominant velocity in x and comparably 
small velocity in y, saying m > with v = £ > 0. The simple arithmetic averaging CT algorithm gives 

^z,i+l/2J+\/2,k - 4 y^:J+l/2.J,k ^:j+\/2J+\.k ^zJJ+l/2,k ^z,,-+i j+i/2,iJ ■ ^-^^f 

* /7+ 1 /2 

In the limiting case of m > with v = £ — > 0, it is obvious that '. ^.^j^j k the only electric field that is in 
the upwind direction, whereas the rest are not. One can easily see that the similar situation also occurs in the 
CT scheme in [30]. This suggests that the simple averaging, based on the ID Riemann solver strategies used 
with CTU, is potentially exposed to numerical oscillations and therefore its stability is questionable. There 
are several other modem time-evolution strategies for MHD that do not suffer from this lack of upwinding. 
As noted, the modem MHD schemes by Balsara [10, 11] using the genuinely multidimensional Riemann 
solvers evolve the magnetic field stmctures in any direction without resorting to any added dissipation in 
the electric fields. The use of multidimensional Riemann solvers for MHD have shown to possess superior 
capability in evolving magnetic fields to the use of conventional ID Riemann solvers, better reflecting the 
tme nature of the PDE that does not require any secondary dissipation mechanisms for the purpose of stable 
upwinding. Although the essential role of the multidimensional technology is acknowledged, our primary 
goal in this paper is to design an easy altemative that can be ameliorated within the ID Riemann solver 
framework based on CTU. 

We now describe our new upwind CT construction scheme that resolves this lack of upwind information 
in the current strategy. As suggested, the idea is to construct the electric fields at {i + + j,k) including 
the electric fields at cell interfaces that are in the upwind directions. For example, in the limiting case of 
u> with V = the cell-comered electric field is given by 

pn+l/2 _ p*,n+l/2 

Z,;+l/2,,/+l/2,i- ~ j,/,,/ + l/2,i- ^■^-^> 

Based on this simple idea of upwinding, we illustrate a systematic approach to constructing a new upwind- 
MEC algorithm that also leads to a consistent CT scheme for plane-parallel, grid-aligned flows. To make 
our discussion more concise, we display a 2D case; the extension in 3D is straightforward. The first step is 
to check the upwind direction at each cell corner. This can be done by defining four switches 

Mp = ^(1 +sign(M/+i/2j+i/2))|sign(M,+i/2j+i/2)|, (56) 

UN = ^(1 -sign(M;+i/2j+i/2))|sign(M,+i/2j+i/2)|, (57) 

vp = ^(1 +sign(v,+i/2j+i/2))|sign(v,+i/2j+i/2)|, (58) 

Viv = ^(1 -sign(v,-+i/2j+i/2))|sign(v,-+i/2j+i/2)|, (59) 
where the sign function is defined by 
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sign(;c) 



1 if .X > 0, 
if x = 0, 
. -1 if ;c<0. 



(60) 



The cell-centered n time step velocity fields are spatially averaged to get the velocities at the cell comer 
(j+ + i) in (56)-(59). When deciding the proper upwind direction at {i + + j) in Equations (56)- 
(59), it is useful to measure relative magnitudes of velocity fields in order to avoid any numerical noise 
effects. The small noise perturbations in the signs of velocity fields may lead to an unnecessary amount 
of changes in upwind directions that are not very advantageous [54]. This motivates to set the velocities 
at (/ + ^,7 + ^) in (56)-(59) to be zero whenever a given local velocity in one direction is relatively small 
compared to the total local velocity magnitude. That is to say, 



M/+l/2,,/+l/2 



if 



/2l 



max 



<ei. 



(61) 



'/+l/2,,/+l/2 ^■+l/2j-+l/2 



Notice that the total velocity only includes the two velocity field components u and v (but not w) that define 
the electric field under consideration. Our choice of an empirically derived value Ei forces to ignore 
any velocity fluctuations that are smaller than 10% of the total magnitude of velocity fields, and set those 
velocities to be zero for determining the proper upwind direction. An arbitrary small value is chosen for £2 
to prevent division by zero. 
Finally we take an upwind biased averaging of the electric fields using the switches in (56)-(59), 



'-'z,i+l/2J+l/2,k 



a 



Vp 



Up 



p*,«+l/2 
^z,i+l/2J,k^ 2 



p,p*,n+l/2 
^y"^z,i+l/2,i.k 



dy 



^y^^zJ+l/2J+l,k , Ay^ ^^^z^i+l/2,j+l,k 



*,n+\/2 



*,n+l/2 

Zj+l/2,j+lk- Y 3^ 
3^*,«+l/2 



. 1 ;^2f*'"+V2 
A/ ^;,/+l/2,y,/ 

8 dy^ 

2 32£*'"+'/2 



+ 



p,*,«+i/2 Ax 

^zJJ+l/21+Y 



z,Li+\/2,k ^ Ax^^^Kl 



dy^ 

*,«+l/2 \ 
j+i/2.k 



dx 



dx^ 



+ 



un E. 



*,n+\/2 
z,i+lj+l/2.k' 



Ax"^z,i+l,j+l/2,k _^ Ax- 



■^axs:£i/2,; 



dx 



3x2 



(62) 



Here the averaging weight factor a is set to 1 if M,+i/2./+i/2i^(+i/2,/+i/2 = 0; a = ^ otherwise. This is our 
upwind-MEC scheme. It is interesting to observe that the upwind-MEC scheme satisfies the consistency 
relationship that reverts to the underlying integration CT scheme for plane-parallel, grid-aligned flows in an 
upwind sense. To see this we consider for example d/dy = 0. Consider first when v,+i/2,/+i/2 = 0. In this 
case the electric fields at each cell corner will take only either the third (if m,+i/2,;+i/2 > 0) or the fourth (if 



M,+i/2,;+i/2 < 0) part of Equation (62). By planar symmetry, E"^ 
terms in both relationships in (62) become 



-1/2 _ 

■j+l/2,k " 



-,H+l/2 _ E,n+l/2 



^z,ij+l,k 



, the first leading 



pn+1/2 

^Z,;+l/2,j+l/2,«: 



^*,n+l/2 _ 
'-'z,i,j+l/2.k ~ 
^*,n+l/2 
'-'z,i+lj+l/2,k 



h,Lj.k if "/+1/2J+1/2 > 0, 

= Ez,i+lJ,k if M,-+i/2j+l/2 < 0. 



(63) 



Therefore they can be considered as an upwind-biased CT variant of 



«+l/2 



■,*,H+l/2 



z,i+l/2J+\/2,k '^zJ+l/2J,k 



which is the 



result of the CT method by Gardiner and Stone [30] in this case. If M/+1/2./+1/2 = then e"^^^^ 



.i+\/2J+\/2.k 
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which is an exact solution for ideal MHD ^ . 
For nonzero values of v,+i/2.,+i/2> consider for instance a case for v/+i/2,/+i/2 > with m,+i/2,/+i/2 > 0- 

Then the electric fields from the upwind-MEC scheme will take the parts that have up and vp only, and 
I 

2- 



a = 5. Consider only the first leading terms in the two parts of Equation (62) for an exposition purpose, we 



get 

pn+\/2 _}_f p*,n+l/2 p*,n+l/2 \ _ }_ f p*,n+l/2 , i7*,n+l/2 \ /^.n 

^z.i+i/2,j+l/2.k — 2 \ z,l+i/2.j,k^^z,i,j+i/2.k j ~ 2 \ z,i+l/2J+l/2,k^^z,iJ+l/2,k j ' 

where we assumed 3/3^ = in the last equality. Compared to the electric field ^"/|'//2 j+i/2 k ~ ^*''i+i/2 j k 
from the Gardiner and Stone's method, the upwind biased MEC electric field in Equation (64) makes use 
of an additional upwind electric field at + ^,k) and includes that field in the average to get the field at 
{i + ^,j + {,k). 

There are several important features of the upwind-MEC method. First, the method appropriately uses an 
upwinding direction rather than taking a simple arithmetic average which lacks proper upwinding. The lack 
of upwinding is found in most of the well known CT schemes [5,30,31,40,47]. The upwinding strategy 
becomes most crucial when advecting a magnetized object in one biased direction, for instance, the weakly 
magnetized field loop advection problem [30] in x-direction only or with a very small advection angle 9 
relative to x-axis. In this small angle advection case the standard CT update without any upwinding becomes 
very vulnerable to numerical instabilities that appear as spurious oscillations in the magnetic field evolution. 
Such oscillations are more likely in the small angle case than in a relatively large angle case because there 
is only one dominating direction from which the CT electric averaging scheme should rely on to obtain 
enough numerical dissipation to stabilize the field evolution. It will be shown later that the upwind-MEC 
strategy advects the field loop without significant numerical oscillations and without distortions for small 
angle advections. 

Second, the upwind-MEC scheme not only accounts for an upwinding direction for stability, but also 
includes high-order terms. The first derivative terms reflect correct spatial changes in expanding from the 
center nodes to the corners, while the second terms effectively avoid spurious oscillations near disconti- 
nuities by adding the proper amount of numerical dissipation to the comer extrapolated fields [40]. These 
high-order terms are upwind averaged in such a way that the scheme is consistent for plane-parallel, grid- 
aligned flows. 

Third, as mentioned, the idea of using upwinding in taking the average is to recover a proper amount of 
numerical dissipation required to ensure stability. We note that the greatest benefit occurs when there is a 
dominating direction locally towards which the magnetic fields are advected. For this reason the upwind- 
MEC scheme can be turned off when the local flow velocities are all ignorable. When the local velocities are 
all negligibly small but finite the local flow should be smooth enough, and hence it is sufficiently accurate to 
use the standard-MFC scheme that takes the arithmetic averaging as discussed in Section 2.6.1. In practice, 
we switch back to the standard-MFC when the local flow velocities are relatively small compared to the 
local sound speed Q. That is, we consider a local Mach number for the local flow switch to choose the 
standard-MFC for constructing the electric field E^^i+i/2.j+i/2,k if 



2 2 
"/+1/2.7+1/2 + ^/+l/2,./+l/2 



M, = ^ ■ — <e3- (65) 



* Note that for non-ideal MHD the upwinding approach in the upwind-MEC scheme should only be applied to the induction part (i.e., — u x B) of 
a generalized Ohm's law including the ternis such as the magnetic diffusion, the Hall effect and the Biermann battery effect. 
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An empirical based tunable parameter £3 = 10""^ suffices to detect a local smooth flow in order to convert 
back to the standard-MEC method; otherwise the upwind-MEC scheme is enabled for all the numerical 
tests presented in this paper. 



2.6.3. CT Update from n to n + I Time Step 

Using the electric fields £^^'^'+^1/2 i:+i/2' ^"/+//2 j k+\/2 ^^'^ ^"t+l^2 j+i/i k constructed by our MEC strat- 
egy, the final CT update evolves the cell face-centered magnetic fields satisfying the V • B = condition on 
a staggered grid. Displaying only in A:-direction, we have 

"x,i+\/2.j.k — "x.i+ll2.j,k 

_AZ_r~„+i/2 _g«+l/2 l_^r_F«+l/2 ,£"+1/2 1 .... 

^yy^z.i+\l2.j+\/2.k ^z.i+l/2.j-l/2,kf ^\ %,/+l/2,;,*:+l/2 + ^y./+l/2,/,A-l/2/- ^""^ 

This completes our description of all procedures in the 3D USM algorithm for a single time step update. 



3. Summary 



We summarize the 3D USM algorithm as follows: 

(i) Calculate the normal predictor states in all .T,3',z-directions using the algorithm described in Section 
2.2. When calculating the normal state in each direction, include the associated MHD term that is pro- 
portional to the gradient of the normal field, see the first relation in (19). During each normal predictor 
calculation, the eigensystems in the normal direction are to be computed. They are stored for later use 
in the transverse correctors. At the same time, the summations of the jumps in all characteristic vari- 
ables are also computed and stored, see Equations (26)-(27), the pseudo-code in Section 2.3.1, the 
sigma summations in Equations (45) and (47). 

(ii) The normal predictor states are updated via the transverse correctors described in Section 2.3.2. This 
step uses two of the stored sigma summation terms that are calculated and stored in each normal 
predictor step. The summations reflect the transverse flux gradients using the characteristic tracing 
approach. 

(a) The reduced 3D CTU scheme then proceeds to advance the normal fields by half a time step 
using CT as illustrated in Section 2.4, finalizing all the interface state calculations. In the reduced 
3D CTU scheme, the formal stability limit is given by a CFL number that is less than j . 

(b) If the full 3D CTU scheme is chosen, the algorithm needs to take one more correction step, 
presented in Section 2.3.3. This correction step is essential in order to provide the full stability 
limit by including the diagonally moving upwind information along the comers in a 3D control 
volume. Similar to the reduced CTU scheme, the full 3D CTU scheme is completed by evolving 
the cell face-centered magnetic fields by half a time step as in Section 2.4. 

So far, both of the two CTU algorithms have required the first set of three Riemann problems that are 
used to advance the magnetic fields by CT. 

(iii) Solve the final set of three Riemann problems at cell interfaces and update the cell-centered conser- 
vative variables to the next time step as described in Section 2.5. The total number of Riemann solves 
therefore becomes six . 

(iv) Calculate the electric fields at cell corners by using the upwind-MEC algorithm described in Section 
2.6.2. With these electric fields, the magnetic fields at cell face centers are updated to the next time 



" Our unsplit data reconstruction-evolution algoritlim can be easily modified for use as a gas liydrodynamics solver by omitting those steps related 
to the magnetic fields. In this case, there are only three Riemann solves required as there is no CT update needed. This unsplit hydrodynamics solver 
has been also available in FLASH'S official releases. 
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step by CT. The cell-centered magnetic fields are updated by taking arithmetic averages of these 
divergence-free magnetic fields at cell face centers (e.g., Equation (21) in Section 2.2). 



4. Numerical Results 



In this section we exhibit the accuracy, stability, convergence and computational performance of the USM 
scheme on a suite of 3D MHD problems. These results show that the scheme is very robust with the full 
CFL stability bound. The scheme is second-order accurate for smooth flows and maintains the solenoidal 
constraint on the magnetic field up to machine round-off error. The full 3D CTU method is our primary 
default method for which a CFL number of 0.95 is chosen in all of the simulations presented here. We 
also show a set of comparison studies between the reduced and full 3D CTU schemes. Insofar as choice of 
Riemann solver is concerned, we use the Roe-type linearized solver [53,63] and the HLLD solver [48]. Our 
choices for the normal predictor step are MUSCL-Hancock, PPM, and WEN05. 



4. 1 . Field Loop Advection 



This problem is notoriously difficult to solve, not because of any strong shock causing numerical instabil- 
ity and leading to code to crash, but rather because it requires full accounting of multidimensional advection 
in a stable matter, such as including the multidimensional MHD terms [30, 31,40, 47]. Failure to do so re- 
sults in an erroneous generation of in-plane magnetic field, which results in the distortion of the initially 
circular (2D) or cylindrical (3D) field loop. 

In addition to the standard field loop advection case studied in [31], we also consider a small angle 
advection case. This turns out to be a much more stringent test than the standard advection configuration 
which assumes a (relatively) large angle between the advection flow and the Euclidean coordinate axes. 
In the small angle configuration there is one dominating coordinate direction along which the field loop 
is advected. This means that, in practice, there is only one direction from which a numerical scheme can 
obtain the numerical dissipation required for stability. In multidimensional problems, inadequate numerical 
dissipation from transverse directions can give rise to anomalous oscillatory behavior in physical variables. 

We begin by describing the initial setup for the standard large angle advection case following the config- 
uration of [31]. The weakly magnetized 3D cylindrical field loop is initialized with a very high plasma beta 
P = p/Bp = 2 X 10^ in the inner region, where Bp = {B\ + B2 +^3)72. Inside the loop the magnetic field 
strength is very weak and the flow dynamics are dominated by the gas pressure. 

The initial field loop is tilted around the X2 (or y) axis by (O = tan~'n radians in a 3D periodic box 
[—0.5,0.5] X [—0.5,0.5] X [—1, 1]. For the standard large angle setup, we choose Q. = 2. The field loop is 
frozen into the ambient plasma and is advected diagonally across the domain with the plasma advection 
velocity (m,v,w) = (1,1,2). The density and pressure are equal to unity everywhere, and Y= f ■ 

The magnetic field components are initialized by taking numerical curl of the magnetic vector potential 
A = (Ai ,A2,A-iY in order to ensure V • B = initially. The relationship between magnetic field and vector 
potential gives 

3A3 3A2 „ 3A3 3Ai „ 3A2 3Ai 
0x2 ax^ 0x1 ox^ ax I 0x2 

For the components of A we choose Aj = A2 = and 

Ao(7?-r) if r<R, 
otherwise. 



(68) 



By using this initialization process divergence-free magnetic fields are well constructed numerically on a 
staggered grid. The parameters in Equation (68) are Aq = 10"-^ and R = 0.3. 
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The two coordinate systems {xi,X2,xj) and {x,y,z) are related by a rotation about the 3'-axis, which is 
given by 



\X3/ 



I cosco sinco ^ 

1 
\^ —sinco cosco j 



(x\ 



y 

VJ 



(69) 



Before we present the standard large angle field loop advection, we first consider a small angle advection 
case in 2D. The 2D initial conditions can be found in [30,40] and we will not repeat the details here. In the 
2D setup, the velocity field is given by 



U = (i^oCOs6,Mosin6, 1)' 



(70) 



where uq = V5. We chose = tan"' (0.01) ^ 0.573° for a small angle advection test. With this setup, the 
field loop is advected almost entirely in the positive .r-direction. This situation makes it hard to stabilize 
the solution during advection because there is not enough numerical dissipation from the j-direction. The 
solution behavior completely relies on the dissipation mechanism from the ;c-direction only, making the 
problem very unstable if no special care is taken to stabilize it. 

In Figure 2, we illustrate the evolution of the magnetic pressure Bp at f = 0.1 and 2.0 using the upwind- 
MEC scheme described in Section 2.6.2. In consequence of the upwinding dissipation mechanism, the 
upwind-MEC scheme stabilizes the solutions extremely well, suppressing the anomalous behavior during 
the advection. 





(a) Bp at ( = 0.1 



(b) Bp at / = 2 



Fig. 2. The 2D field loop advection using a small advection angle 9 0.573° relative to the .v-axis. The images are magnetic pressures at times 
t = 0.1 and 2 using PPM and the Roe Riemann solver. The minmod slope limiter is used for taking slope gradients of characteristic variable in the 
PPM reconstruction step. All results are resolved on 200 x 100 grid cells using the upwind-MEC scheme. 



In Figure 3 (a), the small angle advection test is repeated in 3D, whereas a large angle advection is 
demonstrated in Figure 3 (b). The velocity fields are respectively given by U = (cose, sine, 2)^ and U = 
(1, 1,2)^ for the small and large angle runs. In (a), the same small advection angle ^ 0.573° was used 
relative to the x-axis as in the 2D case. In the large angle case in (b), the field loop makes a domain diagonal 
advection from the given initial velocity condition. We set the tilt angle CO in Equation (69) to be same as 
6 for both (a) and (b). In both runs, we show that the upwind-MEC scheme manifests an oscillation-free 
advections, well-preserving the initial cylindrical shape in the magnetic pressure as shown. As noted above, 
numerical dissipation in the large angle run is naturally added from all directions, rather than from one 
specific direction along the advection in the small angle case. Such added dissipation makes the large angle 
case easier to demonstrate than the small angle case. 
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(a) Bp at/ = 2. (b)Bpat« = 2. (c)B,, atr=l. 

Fig. 3. (a) 3D field loop advection using a small advection angle f» 0.573° relative to the .r-axis. (b) 3D field loop advection using a large advection 
angle using U = (1, 1,2)^. (c) The standard field loop advection problem at time t = I. All results use PPM and the Roe Riemann solver, and the 
minmod slope limiter on characteristic variables in the PPM reconstruction. The results in (a) and (b) are resolved on 64 x 64 x 128 giid cells, and 
(c) on 128 X 128 x 256. The upwind-MEC is used in all cases. 



As a final test in this section, we perfomi the standard field loop advection problem in Figure 3 (c), 
following the configuration in [31]. We see that the upwind-MEC scheme performs very well in evolving 
the field loop successfully to the final time t = I. This result in (c) can be directly compared to the results 
reported in [3 1]. We also report that the upwind-MEC scheme increases the maximum value of the magnetic 
pressure by 48% to 7.41 x 10~^ from its initial value of 5 x 10~^. The larger growth of the maximum value 
is found in the standard-MEC scheme, increasing the initial value by 69% (not shown here). 

Furthermore, we present two quantitative results in Figure 4. They include (a) the temporal evolution of 
the volume-averaged magnetic energy density normalized to the initial (analytic) value < Bp >=< >= 
Bqa/57i7?^/2; and (b) the temporal evolution of the normalized error < \Bi\ > /Bq. Both results in (a) and 
(b) are similar to those reported in [3 1 , 47] . However in (b), the final values at f = 1 are found out to be little 
larger than those in [31, 47] at each grid resolution. This is probably because our full 3D CTU method of 
including the MHD multidimensional terms ignores the 0{At) terms in evaluating the eigensystems of the 
Ad matrices at « + | (see for example Equation (35)), where d = x,y,z- 



4.2. Circularly Pola rized A Ifi'en Wa ve 



In the next test we solve the circularly polarized Alfven wave and its propagation [30, 31,65]. This prob- 
lem provides an important quantitative test of the 3D USM scheme because the smooth initial conditions 
are nonlinear solutions to the problem. The Alfven wave propagates parallel to the xi-axis of a transformed 
coordinate system {xi,X2,xj,) in the periodic computational domain [0,3] x [0, 1.5] x [0, 1.5]. The computa- 
tional domain is resolved on 2N xN xN grid cells, where we adopt = 8, 16, 32 and 64 for the convergence 
study. 

The relationship between the rotated coordinated system {x,y,z) and the non-rotated system {xi,X2,X3) is 
described by the following coordinate transformation 



^ xi\ ( .xcosacosP+j'cosasinP + zsina ^ 



Xl 

\x,J 



— .xsinP + jcos(3 
\^ — xsinacosP — ysinasinP + zcosa J 



(71) 
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Normalized <B„> Evolution 




Normalized <lB3l>/Bo Evolution 



O.IOOOF 



0.0100 r 



0.0010 



0.0001 



0.0 




(a) Normalized, volume averaged magnetic energy density in time 



(b) Normalized 63 error in time 



Fig. 4. Time evolution of (a) the normalized, volume averaged magnetic energy density <Bp >=<B- > and (b) the normalized error < IB3 1 > /Bo- 
Three dilferent results on the grid resolutions of N = 32,64 and 128 are plotted. The full CTU scheme is adopted with CFL=0.95 using PPM and 
the Roe Riemann solver. 



where sin a : 



cosa = and cos(3 — ^ 



The initial conditions we use are the same as the equivalent test problems described in [31]. The initial 
magnetic field is given by 



B = {B,^,B^,,B^y = (l,0.1sin(27ixi/X),0.1cos(27ixi/X))^ 



(72) 



and similarly the velocity field is 



(0,0.1 sin(27lxi A) , 0. 1 cos(27lxi /X))' 
(1,0.1 sin(27lxi /X) ,0.1 cos(27lxi /X) ) ' 



for traveling wave, 
for standing wave. 



(73) 



We set the wavelength X = 1. The density and the gas pressure are initialized by p = 1 and p = 0.1. We 
choose PPM and the HLLD Riemann solver, with the monotonized central (MC) limiter. 

Figure 5 (a) and (b) show the numerical errors on a logarithmic scale obtained with four different grid 
resolutions of = 8, 16, 32 and 64. We test the reduced 3D CTU scheme using CFL=0.475 and the full 3D 
CTU scheme using CFL=0.475 and 0.95 for the convergence study. The errors of the standing and traveling 
waves are plotted in (a) and (b) respectively. For all cases we follow the error calculation formula used by 
Gardiner and Stone [3 1] in order to compare our results with theirs. The results in Figure 5 (a) and (b) show 
a second-order convergence rate of both reduced and full 3D CTU schemes for the smooth Alfven wave 
problem. 

We also measure the relative CPU cost of the full CTU scheme to the reduced CTU scheme, CPUj.gj = 
CPUf_^jy/CPUi--ctu- fii^d that CPUj-gj is about 0.8 on average, which indicates that our full CTU scheme 
with a higher CFL number (e.g., 0.95) is 20% more computationally efficient than the reduced CTU with a 
lower CFL number (e.g., 0.475). The equivalent performance comparison is different in the 6-solve and the 
12-solve algorithms in [31] in that their relative CTU performance turns out to be 1. 
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As the error magnitudes are nearly identical for both standing and traveling wave modes in the reduced 
CTU with CFL=0.475 and the full CTU with CFL=0.95, we conclude that our full 3D CTU scheme exhibits 
better performance while providing numerical solutions that are second-order accurate. 

The figures exhibit a dependence of the truncation error in the full CTU scheme on CFL, in both wave 
modes. The error corresponding to CFL=0.475 is smaller in the standing wave, whereas it is larger in the 
traveling wave simulation. This type of CFL dependence is also seen in [31]. 

Convergence Rate for Standing Wave Convergence Rote for Traveling Wove 




1 10 100 1 1 1 00 

Grid Size Grid Size 

(a) Convergence rate for the standing wave solutions at ^ = 1 .0 (b) Convergence rate for the traveling wave solutions at f = 1 .0 

Fig. 5. The circularly polarized Alfven wave convergence rate for both the standing and traveling wave problems. PPM is used along with the HLLD 
Riemann solver. 



4.3. Orszag-Tang Problem 

The third test problem is the Orszag-Tang MHD vortex problem [49]. We follow the 3D extension [36] 
of the 2D problem in which the initial velocity field is slightly perturbed by £ in the vertical direction. That 
is, the initial velocity field defined on a periodical computational domain [0, 1] x [0, 1] x [0, 1] is written as 

U = (— (1 +£sin27lz)sin27iy, (1 +£sin27iz)sin27i.x,£sin27iz)'^, (74) 

where we use £ = 0.2 as in [36]. The rest are initialized similar to the 2D case so that 

p =7^, = Y, B = (-sin27l3;,sin47U,0)^, (75) 

where y = |. As in the 2D case, the plots in Figure 6 exhibit nonlinear steepening that builds strong discon- 
tinuities from the smooth initial conditions. We show the evolutions of density at f = 0.5 and 1.0 on 128"' 
grid cells. The density images at the top of the domain are very similar to those in the standard 2D case at 
each corresponding time (e.g., see [40]). The flow symmetries are also well preserved in (b) where density 
has developed into more complicated discontinuous flows. The Roe Riemann solver is used with the PPM 
scheme for data reconstruction-evolution in normal direction with MC limiter 
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Max: 6,542 
Mn: 09543 




Z-Axis 



Z-Axis 



Y-Axis 



Y-Axis 



X-Axis 




X-Axis 



(a) Density at t = 0.5 (b) Density at / = 1.0 

Fig. 6. Density plots of tlie Orszag-Tang problem at a resolution of 128^. 

4.4. Rotor Problem 

We extend the 2D rotor problem [5,40, 65] to a 3D case by applying a small velocity perturbation anal- 
ogous to that introduced in the 3D Orszag-Tang problem in Section 4.3. A dense rotating cylinder is ini- 
tialized on a unit cube domain [0, 1] x [0,1] x [0, 1] with non-reflecting boundary conditions. The initial 
velocity field is defined by 



U = (M2rf(l +£sin27lz),V2(/(l +esin27iz),esin27iz)''', 
where 8 = 0.3 and 

-/(r)Mo(3'-0.5)/ro for r < ro 
U2d = -f{r)uo{y-0.5)/r for ro < r < n , 
for r > ri 

'/('")"o(-^-0.5)/ro forr<ro 
V2d = S f{r)uo{x - 0.5) /r for ro < r < n . 
for r > ri 



(76) 



(77) 



(78) 



The density, pressure, magnetic field, and the parameters are initialized as in the standard 2D case given by 

(79) 



' 10 for r < ro 

l+9/(r) forro<r<ri, 
1 for r > ri 



p = l, B = (5/^471,0,0)^, 



(80) 



where uq = 2,ro = 0.1, ri = 0.115,r = ^/{x — 0.5)^ + {y — 0.5)^, and the taper function /(r) is defined by 
/(r) = (ri — r)/ (ri — ro). The value y = 1.4 is used. 

Panels in Figure 7 exhibit contour plots in x-y plane of the density, magnetic pressure and Mach number at 
the final time t = 0. 15. Contour slices are taken at z = 0.5. The problem is solved on a 128^ grid resolution 
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using the Roe solver with PPM. MC limiter is used for the PPM reconstruction. For all cases 40 equally 
spaced contour lines are plotted. All of the contour plots show that our 3D results correspond very closely 
to the underlying 2D solutions (e.g., see [40]). As reported in [65] one important feature to observe in this 
problem is to check the oval contours of Mach number near the center. As illustrated, the contour lines are 
symmetrical and well preserved with our choice of CFL=0.95. 




0.2 0.4 0.6 O.B 0.2 0.4 0.6 .B 

X-Axis X-Axis 



(a) Density contour at ; = 0.5 ranging between 0.4540 and 14.82 (b) Magnetic pressure contour at z = 0.5 ranging between 0.009705 and 3.171 




0.2 0.4 0.6 O.B 

X-Axis 



(c) Mach number contour at ; = 0.5 ranging between 3.268 x 10"^^ and 5.938 
Fig. 7. The rotor problem with a resolution of 128' at f = 0.15. In all cases, 40 equally spaced contour lines are plotted. 
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4.5. Cloud & Shock Interaction 



In the next test problem we consider the interaction of a high density cloud with a strong shock wave, 
originally studied by Dai and Woodward [19] and often referred to as the cloud-shock interaction problem. 
This problem aims to test the code robustness in solving flow conditions such as high supersonic Mach 
numbers in the pre-shock and the post-shock regions, wide ranges of plasma beta values across the front/rear 
of the cloud, and strong shear flows [36, 40,47, 65]. 

Our computational domain is a cube, spanning from -0.5 to 0.5 in all three directions and is resolved on 
128-^ grid cells. Supersonic inflow boundary conditions are imposed at the lower boundary x = —0.5, while 
outflow conditions are used elsewhere. The initial condition has different left and right states, separated by 
an initial discontinuity at .x = 0. 1, given by 

^ _ /(3-86859,0,0,0,0,2. 1826182,-2.1826182, 167.345) if;c<0.1, 
{p,u,v,w,B,,By,B,,p) - I 2536,0,0,0,0.56418958,0.56418958, 1) if x > 0.1. ^ ^ 

The high density cloud is located on the right side of the domain, and has a spherical envelope defined by 
{x — 0.3)^ + 3^- = 0.15^. A uniform density p = 10 and pressure p = 1 are fixed in the inner region 
of the cloud, and y = 5/3 is used everywhere. The velocity and the magnetic fields are the same as the 
surrounding right state plasma values. The simulation is carried out to a final time t = 0.06 using the 
WEN05 reconstruction scheme and the Roe Riemann solver. We used the van Leer's slope limiter for 
limiting characteristic variables in the WEN05 reconstruction. 

Figure 8 shows density (plotted in the top half using a red color scheme) and magnetic pressure (plotted 
in the bottom half using a blue color scheme) att = 0.06. The main features of the cloud-shock interaction 
process are well captured, in that the temporal evolution of the high density cloud produces disrupted shapes 
as the cloud moves into the plane shock on the left. 

Simulations of this problem are often performed using a rather diffusive set of numerical options such as 
the minmod slope limiter, HLL-type Riemann solvers, or lower values of CFL. For example, as noted by 
Toth, dimensionally-split MHD algorithm can easily fail due to unphysical states (e.g., negative pressure 
or density) arising in consequence of the strong interaction between the shock and the cloud. By contrast, 
the 3D USM scheme utilizing the full 3D CTU algorithm and the upwind-MEC scheme can run this simu- 
lation successfully without relying on such numerically dissipative choices. Despite our choice of the van 
Leer's limiter for WEN05, and of the Roe Riemann solver using CFL=0.95, the final time step is reached 
successfully without giving rise to any numerical instabilities. 



4.6. MHD Blast Wave 



The last test case is the 2D MHD spherical blast wave problem of Zachary et al. [68]. We presented our 
2D results in [40] and extend the problem to 3D here. We test three different configurations, differing by 
the initial strength of magnetic field in x-direction, each leading to strong shock formation and propagation. 

The computational domain is a unit cube [—0.5,0.5] x [—0.5,0.5] x [—0.5,0.5] with a grid resolution of 
128^. The ambient gas is initialized as 

P = 1,P = 0.1,B = (B,„,0,0)^ (82) 

where the three simulations have initial values of B^^ given by B^o = 0> ^.vo = and Bv„ = At the 
center of the domain, a spherical region of radius r = 0. 1 is initialized with a very strong pressure p = 1000. 
The non-zero values of Bx„ = and produce very low-P ambient plasma states, |3 = 1 x 10~^ and 
2.5 13 X lO"'* respectively. Through these low-|3 ambient states, the explosion initially emits almost spherical 
fast magneto-sonic shocks that propagate with the fastest wave speed. The flow has y = 1 .4. 
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(a) Density and magnetic pressure at / = 0.06 



Fig. 8. Tlie 3D MHD interaction between the liigli density cloud and shock structures resolved on 128^ grid using the Roe Riemann solver and the 
5th order WENO scheme. Plotted are density (denoted as "dens" in the legend) in the top half and magnetic pressure (denoted as "magp") in the 
bottom half. 

Shown in Figures 9-11 are (a) density (plotted in the top half) and magnetic pressure (plotted in the 
bottom half) and (b) contour plots of gas pressure (top half) and total velocity U = Vu^ + + (bottom 
half) at time t = 0.01. The contour slice plots are the x-y planes taken at z = 0. 

This problem is susceptible to a type of shock wave instability known as the carbuncle phenomenon [52]. 
The carbuncle instability takes place in multidimensional numerical solutions when using a less dissipative, 
ID based (rather than the multidimensional based [10]) Roe-type Riemann solver, in the regions where a 
planar shock is aligned to the grid. The cause of this instability is the lack of numerical diffusivity added 
to the Roe-type fluxes perpendicular to the grid-aligned shock, resulting in a growth of small amplitude 
noise in the transverse direction. There are several approaches to fix the instability [33,50,56,59] which 
all basically provide a similar mechanism to add extra numerical diffusion in the transverse direction. Here 
we use a hybrid Riemann solver that appropriately combines Roe and HLLE depending on the strength 
of shocks. In this approach, the HLLE solver is adaptively used only in strong shock fronts detected by a 
shock switch [5]; the Roe solver is used elsewhere. The second-order accurate MUSCL-Hancock scheme 
is used for the normal predictor calculations. We also employ a hybrid-type of slope limiter that combines 
MC limiter for linearly degenerate waves (i.e., Alfven and entropy waves) and the minmod limiter for 
genuinely nonlinear waves (i.e., magneto-sonic fast and slow waves). This hybrid limiter approach provides 
an added robustness and accuracy by using a compressive limiter (such as MC and van Leer's) for crisper 
representation of the linear waves, whereas a diffusive limiter (such as minmod) for the self-steepening 
nonlinear waves [7]. 

The case B^^ = is illustrated in Figure 9. The carbuncle phenomenon can appear to be stronger in this 
hydrodynamic limit than when B^^^ ^ 0. Using the hybrid Riemann solver, however, we do not see any 
artifacts at the shock fronts that are aligned to the grid axes. In the absence of magnetic field the explosion 
propagates the shock wave spherically in all radial directions, as exhibited in the contour plots in Figure 9 



In Figure 10 the intermediate magnetic field strength case with B^ = is shown. The explosion becomes 



(b). 
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(a) Density and magnetic pressure at / = 0.01 



(b) Contours of gas pressure and total velocity at f = 0.01 in the x-y plane at 

.' = 



Fig. 9. Results of the blast problem simulation with By = using a hybiid Riemann solver. In (a), density (denoted as "dens" in the legend) is 
plotted at the top half. Magnetic pressure (denoted as "magp" in the legend) is plotted at the bottom half and is represented as small values that are 
numerical noise. In (b), 40 contour lines are plotted for gas pressure (top half) between 0. 1 and 73.62 and total velocity (bottom half) between and 
8.810. 

anisotropic because of the non-zero magnetic field strength in .x-direction. The intermediate value of Bjc 
still permits shock wave propagation in the j-direction, so that the overall spherical shape is not radically 
distorted. Nonetheless, the development of the elongated wave structures in the direction parallel to the 
field is evident compared to the hydrodynamic limit case in Figure 9. 

Finally, Figure 11 illustrates the strongest magnetic field case, = The explosion now becomes 
highly anisotropic. This strong anisotropic behavior is well shown in Figure 1 1(b) in that the displacement 
of gas in the transverse j-direction is increasingly inhibited and hydrodynamical shocks propagate almost 
entirely in the ;c-direction parallel to Bj^. It is also evident that several weak magneto-sonic waves are radiated 
transverse to the A:-direction. This process continues until total pressure equilibrium is reached in the central 
region. 

Balsara [7] pointed out that maintaining positivity of pressure becomes challenging due to the strong 
wave propagation oblique to the mesh. Such unphysical pressures can distort contours, especially near the 
outer boundary where a large and unphysical drop in pressure takes place immediately ahead of the shock. 
In our calculations, the pressure remains always positive throughout the simulations, evidence that our 3D 
MHD scheme is very robust and accurate with our choice of high CFL=0.95. 

5. Conclusion 

We summarize several key features described in this paper. First, the 3D USM scheme has been intro- 
duced, developed and studied. The method is a 3D extension of the 2D USM algorithm [40] which employs 
characteristic analysis to account for contributions of both normal and transverse MHD fluxes in a truly 
unsplit fashion. Therefore they do not need intermediate Riemann solves to correct the normal predictor 
states with the transverse flux updates as in the usual 12-solve CTU algorithm [47,55]. Our approach of 
using characteristic analysis provides computational efficiency by storing the eigensystem evaluations when 
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(a) Density and magnetic pressure at f = 0.01 (b) Contours of gas pressure and total velocity at f = 0.01 in the x-y plane at 

z = 

Fig. 10. Results of the blast problem simulation with 6, = using a hybrid Riemann solver. In (a), density (denoted as "dens" in the legend) is 
plotted at the top half. Magnetic pressure (denoted as "magp" in the legend) is plotted at the bottom half. In (b), 40 contour Hnes are plotted for gas 
pressure (top half) between 0.009981 and 106.6 and total velocity (bottom half) between and 11.84. 




(a) Density and magnetic pressure at f = 0.01 (b) Contours of gas pressure and total velocity at f = 0.01 in the x-y plane at 

z = Q 

Fig. 1 1 . Results of the blast problem simulation with 6, = using a hybrid Riemann solver. In (a), density (denoted as "dens" in the legend) is 

\/47I 

plotted at the top half. Magnetic pressure (denoted as "magp" in the legend) is plotted at the bottom half. In (b), 40 contour lines are plotted for gas 
pressure (top half) between 0.009161 and 202.9 and total velocity (bottom half) between and 15.97. 
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computing each nomial direction and re-using tliem in tlie transverse flux calculations. 

We introduced two different methods, the reduced and full 3D CTU schemes. The reduced scheme can 
be considered as a straightforward extension of the 2D CTU algorithm of Colella [16], and is analogous to 
the 6-solve algorithm in 3D by Gardiner and Stone [31]. Although the reduced CTU scheme has a simple 
implementation for 3D, its stability limit is bounded by a CFL number less than 0.5. The full CTU scheme 
significantly improves this limited stability range and can utilize the maximum stability range of CFL num- 
ber close to 1 . This was achieved by taking into account the second- and third-order cross derivative terms in 
computing intermediate states at n + | and n + j. Our full CTU scheme thus includes the multidimensional 
upwind information that is crucial to provide the full CFL hmit. We also showed that the relative CPU cost 
of the full scheme compared to the reduced scheme is less than 1, indicating the cost efficiency of the full 
CTU scheme. The multidimensional MHD terms are also properly included in both normal and transverse 
directions. 

Second, we extensively investigated the lack of numerical dissipation mechanisms in the existing CT 
algorithms, especially when there is a biased direction in advecting magnetic fields. In the small angle 
advection tests in 2D and 3D, we showed that the field loop simply can fail to be cleanly advected, and 
become distorted into non-circular or non-cylindrical shapes in most CT schemes. By contrast, the upwind- 
MEC scheme, by incorporating upwind information adds the needed numerical dissipation when taking the 
arithmetic average in CT. The algorithm enhances the previous MFC scheme [40] in that upwind-MFC 
maintains consistency of plane-parallel and grid-aligned flows [30]. 

The results of the test problems in Section 4 give considerable confidence in our scheme for use as a robust 
and reliable second-order, finite-volume 3D MHD algorithm. The methods developed in this paper for the 
3D USM scheme preserve the divergence-free constraint without any evidence of numerical instability or 
accumulation of unphysical errors using a very high CFL number close to 1. The suite of test problems 
presented in this study include several stringent setups can be particularly challenging for MHD algorithms. 
The scheme has been thoroughly tested and has been shown to perform very well, providing confidence in 
its ability to correctly simulate a wide range of MHD phenomena. 

The 3D USM scheme presented here has been implemented on both uniform and AMR grids. It has been 
integrated and tested in the official FLASH4 release of the Flash Center for Computational Science at the 
University of Chicago [29]. 
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